Photochemical and thermochemical pathways to S2 and polysulfur formation in the atmosphere of Venus

Polysulfur species have been proposed to be the unknown near-UV absorber in the atmosphere of Venus. Recent work argues that photolysis of one of the (SO)2 isomers, cis-OSSO, directly yields S2 with a branching ratio of about 10%. If correct, this pathway dominates polysulfur formation by several orders of magnitude, and by addition reactions yields significant quantities of S3, S4, and S8. We report here the results of high-level ab-initio quantum-chemistry computations that demonstrate that S2 is not a product in cis-OSSO photolysis. Instead, we establish a novel mechanism in which S2 is formed in a two-step process. Firstly, the intermediate S2O is produced by the coupling between the S and Cl atmospheric chemistries (in particular, SO reaction with ClS) and in a lesser extension by O-abstraction reactions from cis-OSSO. Secondly, S2O reacts with SO. This modified chemistry yields S2 and subsequent polysulfur abundances comparable to the photolytic cis-OSSO mechanism through a more plausible pathway. Ab initio quantification of the photodissociations at play fills a critical data void in current atmospheric models of Venus.


Table of Contents
. Photochemical outcome of the cis-OSSO and trans-OSSO runs departing from the S2 state for the MS-CASPT2 NAMD simulations.

S13
Supplementary Table 2. Information of the NAMD simulations run at the TD-B3LYP/6-31G* level of theory.

S14
Supplementary Table 3. TD-DFT dissociation limits for the (SO)2 studied isomers. S15 Supplementary Table 4. Summary of the number of trajectories per channel and the calculated channel yields for the TD-DFT NAMD runs.   Table 6. MS-CASPT2 activation energies (∆E ‡ ), energy difference between reactants and products (∆E), and calculated rates for thermal processes derived from the reactivity between 1 OSSO and 2 NO / 3 O / 3 S / 2 H, as well as those related to the reaction between 2 ClS and 3 SO. Energetic profiles of the association of 3 S and 3 O2 / S are also shown. Table 7. Photolysis rates from cis-/trans-OSSO to produce SO from Frandsen et al., 2 with Venusian altitude set to 64 km and latitude 0º and updated value considering the yields of photogenerated SO computed in this work.

S35
Supplementary Table 8. Number density (reported and computed here; in molecules cm -3 ) in the Venusian atmosphere at 64 km altitude of the species involved in the bimolecular reactions studied in this work. Table 9. Electronic activation energies (∆E ‡ ) and energy differences between reactants and products (∆E) for reactivity derived from the interaction between 1 OSSO and 3 SO.

S51
Supplementary Table 10. Gibbs activation energies (∆G ‡ ) and Gibbs energy difference between reactants and products (∆G) for reactivity derived from the interaction between 1 OSSO and 3 SO.

Supplementary Note Computational details of the MRCI profiles
The excited-state calculations were done in the Cs symmetry group. We have used the completeactive-space self-consistent field (CASSCF) method followed by the internally contracted multireference configuration interaction (MRCI) 4,5 scheme to compute the dynamic electron correlation.
Here, the atoms were described by aug-cc-pV(T+d)Z basis set. The CASSCF active space was chosen after considering the [1a′-14a′] and [1a′′-2a′′] lowest molecular orbitals as doubly occupied and keeping the remaining valence orbitals as active. MRCI calculations included all configurations of the CI expansion of the CASSCF wavefunctions having a weight larger than 0.01. All singlet and triplet electronic states were calculated to the first dissociation limit (SO + SO), i.e. four 1 A′ and two 1 A′′ states were averaged in the CASSCF procedure. Only three triplet states were averaged, i.e. 3 A′ and two 3 A′′ states. This leads for instance to considering more than 10 8 contracted and 3·10 11 uncontracted configuration state function when computing singlet electronic states. The MOLPRO2019 software was used 6 .

Simulation details
All trajectories were computed making use of the surface-hopping including arbitrary couplings (SHARC) code 7,8 . The electronic structure was computed with the MS-CASPT2 method 9-11 as implemented in the OpenMolcas software package 12 , interfaced with the SHARC program. Therefore, wave functions, energies, nuclear and electronic gradients and spin-orbit couplings were calculated on the fly by the OpenMolcas software whereas SHARC computed the diabatic states, the hop probabilities between them, and propagated the excited-state trajectories. Decoherence corrections were included using the energy-based method of Granucci, Persico and Zoccante with the recommended parameter of C = 0.1 a.u. 13 For all systems, the integration of the nuclear motion was done by means of the Velocity-Verlet algorithm using a time step of 0.5 fs.
Initial conditions were obtained stochastically at 300 K by sampling the ground state minimum via a Wigner distribution employing CASPT2 frequencies. The atomic coordinates and velocities were subsequently used to start the molecular dynamic simulations. For cis-OSSO, 30 trajectories departed from the S1 state and 86 from S2, whereas for trans-OSSO 70 runs started from S2. The S1 state is mostly dark and unreactive. The dynamics were run on diagonal potential energy surfaces (PESs) resulting from the diagonalization of the Hamiltonian containing non-adiabatic couplings and spinorbit couplings 14 , employing the SHARC algorithm for surface hopping. Globally, the total simulation time was 140 fs, even though some trajectories were extended to further simulation times to confirm certain photoproducts. A total of 4 singlet and 4 triplet spin-free (SF) states were mixed to obtain the corresponding SO states, making use of the state-average CASSCF method and employing the CAS(12,10) active space shown in Supplementary Figure 1. All states were active in the simulations except the fourth triplet root, which was excluded in all cases to save computational time.
The MS-CASPT2 energies and gradients were computed with the IPEA shift set to 0.25 a.u. 15 , whereas the imaginary level shift was set to 0.2 a.u. to minimize the presence of weakly-interacting intruder states 16 . The ANO-S-VDZP basis set was used for all MS-CASPT2 dynamic runs.
A certain number of trajectories were excluded from the final analyses (100 trajectories were prepared in the initial setup for the simulations starting from S2) due to premature crashes caused mostly by errors in the calculation of the MS-CASPT2 numerical gradients, problems with SOC determinations, divergence in CASSCF iterations, or inconsistencies in the total energy (variations larger than 0.5 eV between consecutive steps).
As previously used in the literature 17,18 , bonds are considered as broken when they are stretched 1.5 times the bond distance at the Franck-Condon geometry (optimized at the same level of theory as the dynamics). These values are indicated with grey solid lines in the distance plots (also for the TD-DFT NAMD simulations, see below). This criterion is again confirmed by the fact that, considering the significant bond stretching momentum observed in the dissociations, as shown in the corresponding interatomic distance curves, and the fact that for all bonds considered in this work, no trajectory reduces the bond distance after passing the 1.

Method validation
The MS-CASPT2/ANO-S-VDZP method detailed in the previous section has been validated by comparing the energy profiles along the rigid scans of the S-S and O-S bonds for both cis-and trans-OSSO, computed at this level of theory and with the MRCI method ( Supplementary Figures 2-5). Overall, MS-CASPT2 reproduces the main features of the S0-S2 and T1-T3 curves with maximum differences of a few tenths of eV. Discrepancies between both electronic structure methods increase for the T1 and T2 curves at large S-S distances, however, the population of triplet states is clearly negligible and thus these differences will not influence the dynamics. In addition, at these regions, the molecule is already dissociated and therefore the mentioned small changes in energy do not affect the branching into the different photoproducts. Globally, these results fully validate the use of the MS-CASPT2/ANO-S-VDZP method as a suitable level of theory to run the NAMD simulations for both OSSO isomers.

Results
For cis-OSSO and trans-OSSO, two excited states were initially populated, namely the dark S1 and the bright S2 states. Results  For cis-OSSO, about 60% of the S2 population decays into S1 before 40 fs, while the same phenomenon takes place in trans-OSSO at about 80 fs. The triplet state participation is almost zero. In both isomers the S-S photolysis takes place in the S2 and S1 states, the ground state only starts to be marginally populated (about 15%) after 80-100 fs. The SO fragments are thus released in the singlet manifold (S2, S1, and S0 states). Those excited 1 SO fragments are assumed to decay finally to their respective triplet ground states 3 SO. To quantify the 1 SO/ 3 SO ratio, photodynamics should be extended far beyond the time limit of the current simulations, which is prohibitive at the level of theory required for a correct description. This is out of the scope of this work. Nevertheless, future implementations of machine learning techniques together with the photodynamics methodology shall help to quantify this ratio and the relevance of further reactivity of 1 SO.  3. Supplementary Note 3. TD-DFT non-adiabatic molecular dynamics

Simulation details
The simulation protocol for the TD-DFT 19,20 non-adiabatic dynamics is analogous to the one used for the MS-CASPT2 simulations, only the differences between them are mentioned here. Electronic energies and analytical gradients are computed at the TD-B3LYP/6-31G* level 21,22 with the Gaussian 09 program 23 , whereas SHARC computed the diabatic states, the hop probabilities between them, and propagated the excited-state trajectories. The ground state was computed using a restricted DFT ansatz, and the quadratic convergence was used when default self-consistent field procedure did not converge. No triplet states were included in the simulations since its participation is excluded from the MS-CASPT2 simulations.
Initial conditions were generated with a Wigner distribution as commented before, however, the population of the initial state was decided through a different procedure. We performed single-point calculations with the TD-B3LYP/6-31G* method on top of each initial condition to obtain a set of vertical absorption energies and oscillator strengths. Later, the initial active state was flagged stochastically, as documented elsewhere 24 , simulating an excitation window that varied according to the optical properties of the different isomers. Supplementary Table 2 collects the excitation windows, number of trajectories, states included in the simulations, initial active states, and total simulation time for the NAMD at the TD-DFT level. The simulations were run only for 80 fs, since it is sufficient time to observe the photodissociations (and the mutual competition between the different breaking channels) and the population of the ground state (not well accounted in the TD-B3LYP NAMD, see next section) is negligible, as shown by the population analysis of the MS-CASPT2 simulations ( Supplementary Figures 8 and 10).

Method validation
This section explains in detail the validation of the NAMD dynamics performed at the TD-B3LYP/6-31G* level. This level of theory has been systematically benchmarked against MS-CASPT2 dynamics and static profiles to ensure the full reliability of the obtained results.
In bonded structures (i.e. when the molecule is not fragmented), results clearly show that the TD-B3LYP/6-31G* method accurately reproduces the excited state energy ordering and associated oscillator strengths as compared to MS-CASPT2(12,10)/ANO-S-VTZP determinations. This is clearly shown for several 1 (SO)2 isomers in Supplementary Figures 11-16. Therefore, there is a range of bond distances, including regions close to the Franck-Condon area and zones with relatively long bond lengths (see below), in which the TD-B3LYP/6-31G* description is fully reliable. Since the branching ratio into different photodissociation channels is largely dominated by the initial coordinates and velocities and the nature of the initial active state, and considering that the ultrafast dissociations take place at the sub-100 fs regime, where the ground state is not active at all, this level S15 of theory provides the correct photoproduct distribution for each system, which is the main goal of this study. As a matter of fact, TD-DFT dynamics correctly identify the 3 SO + 3 SO dissociation as the only photoproduct (100%) of cis-and trans-OSSO, in reasonable agreement with the MS-CASPT2 dynamics (≥90%).
On the other hand, it is well known that the ground-state profiles computed with the restricted DFT ansatz are not reliable at bond dissociation limits because a single restricted determinant cannot account for both zwitterionic and diradical solutions (see top panels of Supplementary Figures 11-16). Thus, these methods typically show an artificial energy increase as long as the bond distance elongates, whereas multiconfigurational schemes such as MRCI or MS-CASPT2 provide the correct asymptotic, smooth dissociation profiles. Thus, the bond distance at which the energy starts to increase artificially, can be easily found by comparing the TD-B3LYP vs MS-CASPT2 profiles. This has been systematically evaluated for all systems considered in this work in Supplementary Figures 11-16. By inspection of these results, we can clearly see that, for cis-OSSO, the TD-DFT profiles exhibit an artificial minimum at a S1-S3 distance of ~2.8 Å followed by a significant energy increase upon larger interatomic distances, whereas the MS-CASPT2 profiles are more planar (Supplementary Figure 11). Therefore, this minimum can artificially trap some trajectories in this area that would have certainly dissociated at the MS-CASPT2 level, as unambiguously observed in the MS-CASPT2 dynamics (Supplementary Figure 7). The reason is that once reached these S1-S3 bond distances, the molecules always dissociate. Therefore, it is reasonable to conclude that the runs that reach this area must be considered as dissociated. Taking into account that the description at shorter distances is trustworthy (as shown by Supplementary Figure 11), it can be safely concluded that the dissociation limit for the S1-S3 bond at the TD-DFT level is 2.8 Å. This value is indicated with a thicker black line in the corresponding figures of the next subsection.
The same analysis can be performed for the rest of systems and bonds. Supplementary Table 3 summarizes the TD-DFT bond dissociation limits determined by comparison with the MS-CASPT2 profiles. Table 3. TD-DFT dissociation limits for the (SO)2 studied isomers.

System
Bond Photoproduct Finally, in order to further test the TD-B3LYP NAMD validity, we have computed the potential energy landscape along a cyclic-OS(=O)S → S + SO2 trajectory with both the TD-B3LYP/6-31G* and the MS-CASPT2(12,10)/ANO-S-VTZP methods. Results are shown in Supplementary Figure  17. The agreement between both profiles is very good, especially at the beginning of the run.
Three important facts emerge from all MS-CASPT2 and TD-DFT analyses presented above: S16 • The MS-CASPT2 dynamics indicate that neither the ground state nor the triplet states participate in the photodissociations of the cis and trans-OSSO. • In TD-B3LYP NAMD, hops to the ground state are almost shut down due to the poor coupling between the ground state (single determinant) and the excited states (CI expansion of single excitations and de-excitations) due to the TD-DFT formalism. Therefore, artificial trapping in ground/excited-state degeneracies (not happening at the MS-CASPT2 level, in which both ground and excited states are treated on equal foot and normal hopping is allowed), are not a problem. • The dissociation bond distances can be easily corrected by comparing both TD-B3LYP/6-31G* vs MS-CASPT2/ANO-S-VTZP energy profiles.
The following reasons justify the use of the TD-B3LYP/6-31G* method to describe the photodynamics of the 1 (SO)2 isomers. Note that other molecules not studied here would require the corresponding validation studies.
• The description of the most relevant excited states and oscillator strengths at the TD-B3LYP/6-31G* method is excellent at distances shorter than the dissociation limits, ensuring the correct branching into the different photodissociation channels and therefore allowing correct quantifications of the photoproducts.

Supplementary Note Static ground-state CASPT2/MS-CASPT2 reactivity 4.1 Computational details and benchmark analyses
The chemical reactions summarized on Table 1 of the main text were modeled by means of CASPT2 optimizations in combination with the ANO-L-VTZP basis set. The reason for the choice of thermal reactions was to improve old rough estimations of kinetic rates involved in the direct or stepwise generation of S2 in an alternative manner as photochemical pathway proposed by Pinto et al. 1 (reaction 2 of the main manuscript). We update in Supplementary Table 6 rates related to reactions with energy barriers using conventional transition state theory. Note that for barrierless association reactions forming (SO)2 dimers, such as 3 SO+ 3 SO, which are also relevant in the overall modeling, Frandsen et al. 2,25 have recently reported computational data based on single-reference methods and collision theory and variational transition state theory. We used the literature rates for the barrierless reactivity.
Minima and transition states were optimized numerically using the subroutines implemented in OpenMolcas 12 . The true nature of the TSs was verified by analyzing the list of vibrational frequencies (only one negative vibration corresponding to the reaction coordinate) and through intrinsic reaction coordinates (IRCs) calculations to connect the TSs with the corresponding reactants and products. These latter structures were further optimized without constraints to ensure the full relaxation of the structures. In some cases, relaxed scan calculations and/or linear interpolation in internal coordinates (LIIC) were also used to approximate the reaction profiles, using an in-house code for the interpolations and again OpenMolcas to compute wave functions and energies on top of the interpolated structures.
The different active spaces used for each ground-state reaction are listed on Supplementary Table  5. All CASPT2 optimizations were performed computing only one root (state-specific) in the CASSCF procedure. To test the influence of averaging more states in the CASSCF method in the reaction profiles, state-average (SA)-CASSCF computations including 3 roots were also performed on top of the converged structures. The dynamic electron correlation of the 3 states was recovered with the MS-CASPT2 ansatz. With this latter protocol, namely SA(3)-MS-CASPT2, only the first root is displayed in the reaction profiles.
Results show that there is significant state-average effect in some systems, especially in regions close to the TS structures in which the ground and excited states couple more efficiently. A comparative analysis between the CASSCF and the SA(3)-CASSCF wave functions reveal slight differences in the weights of some important configuration state functions, causing the energy differences. In general, the energy barrier heights tend to be smaller at the SA(3)-MS-CASPT2 level of theory as compared to the state-specific CASPT2 method. For reaction (3), in which an accurate experimental rate constant is reported in the literature (~7.5 kcal/mol) 26 , the MS-CASPT2 data for reaction (3) gives an energy barrier of 9.0 kcal/mol, showing a clearly better agreement with the experiment. Therefore, SA(3)-MS-CASPT2 energy profiles must be considered to have the highest accuracy. State-specific CASPT2 profiles are used for comparison purposes to display the effects state average and multistate approach.   Fig. 41 a At the intermolecular distances of the post-reactive complex, the doublet or triplet state is delocalized over the two reaction products. Table 6. MS-CASPT2 activation energies (∆E ‡ ) and energy difference between reactants and products (∆E) and calculated rates for thermal processes derived from the reactivity between 1 OSSO and 2 NO / 3 O / 3 S / 2 H, as well as those related to the reaction between 2 ClS and 3 SO. Energetic profiles of the association of 3 S and 3 O2 / S are also shown. The value estimated from experimental measurements for the reaction 3 SO + 1 OSSO → 1 SO2 + 3 S2O is shown within parenthesis. Unimolecular and bimolecular rate constants have been computed by means of transition state theory at T=298 K. Pseudo-first order rates were calculated using the concentration values of the reactants presented in Supplementary

Supplementary Note Excited-state MS-CASPT2 energy profiles
The energetic profiles for some particularly interesting excited-state dissociations have been computed with the MS-CASPT2/ANO-S-VDZP method and the OpenMolcas software. 12 Figure 42. MS-CASPT2 relaxed scan of 1 cis-SSSO, relaxing the S2 state. The absorption energy to populate the S2 state at the S0 min region (S1-S3 distance of 2.16 Å) is shown through the dotted arrow. The energy of the S0 and S1 states at the S1 min geometry (S1-S3 distance of 2.21 Å) is also shown. Four roots have been computed in the SA-CASSCF method. 6. Supplementary Note 6. Ground-state DFT reactivity 6

.1 Computational details
DFT optimizations were carried out using the B3LYP/6-311G(d,p) level of theory as implemented in the Gaussian 16 program. 29 The restricted scheme was employed for systems with singlet multiplicity, whereas the unrestricted DFT ansatz was employed for triplet states. The nature of each stationary point was verified by examining the list of frequencies (all positive for minima, one negative and the rest positive for transition states). Intrinsic reaction coordinate (IRC) computations were systematically computed to ensure the connectivity between reactants, transition states, and products. Unconstrained optimizations were performed to the final structures of the IRC calculations in order to ensure the complete relaxation of reactants and products.

Results
In Supplementary Figure  The formation of TS2 implies the breaking of the S2-S3 and O6-S bonds, giving a structure with a relative energy of -0.03 eV. The energy barrier from I1 to TS2 is thus of 0.46 eV (10.6 kcal/mol). The release of the 1 SO2 molecule and the stabilization of the 1 cis-S3O2 is barrierless, releasing 0.56 eV (12.9 kcal/mol) and yielding the van der Waals complex P1 ( 1 SO2 + 1 cis-S3O2) with a relative energy of -0.59 eV (13.6 kcal/mol).
The left-hand side of Supplementary Figure 45 shows the reaction 1 cis-S3O2 + 1 SO2 → 1 S2 + 2 1 SO2. The 1 cis-S3O2 decomposition to yield the products P2 has a global energy barrier of 0.81 eV (18.7 kcal/mol) mainly associated to the formation of the linear species S-S-O-S-O (I2). The chemical channel goes through two planar areas, TS3 and TS4. The bond breaking of the O8-S7 that releases the second SO2 molecule and S2 is barrierless. Further reaction profiles with SO are shown in Supplementary Figures 46-47. On the other hand, Supplementary Figures 48-50 compute the reaction profiles for an alternative reaction of 1 cis-S3O2 with 3 SO.

S47
Supplementary Figure 45. Reaction of two 1 cis-OSSO molecules to yield 1 cis-S3O2 and 1 SO2 (right), followed by the 1 cis-S3O2 decomposition to yield 1 SO2 and 1 S2 (left). Reactants are at the right-hand and products at the left-hand of the plot in coherence to the progression of the S3-O8 bond shortening. The profiles have been obtained by means of TS optimizations followed by IRC calculations to connect reactants, TSs, and products (right) and through relaxed scans of the S3-O3 bond distance (left) shrinkages at the B3LYP/ 6-311G(d,p) level of theory. Dashed lines indicate non-connected points due to molecular reorganizations involving other molecular coordinates that have been studied separately (see below). TS3 and TS4 areas refer to transition state areas (shadowed in the plot) identified along the PES by visual inspection, the highest-energy structures are not the result of TS optimizations.  Figure 47. Exploration of the TS4 region through a combination of LIIC (right) and MEP (left) techniques. The PES is planar along the S3-O8 bond shortening from 1.644 to 1.544 Å, as shown by the LIIC profile, with a relative energy of -1.34 eV with respect to the energy of the reactants. The MEP indicates that the formation of the second 1 SO2 molecule and the final release of singlet diatomic sulfur 1 S2 is barrierless. This profile leads to the products P2 (2 1 SO2 + 1 S2).

Computational details
The ground-state reaction profiles computed in the previous section were recomputed through optimizations at the CCSD/cc-pVDZ level, employing Gaussian 16 29 . Frequency calculations served to identify the nature of the stationary points. IRC calculations were performed at the same level of theory using the ORCA 4.2 software 30 .
Single-point corrections using the CCSD(T) method were computed in combination with two larger basis sets, namely the Pople 6-311+G(3df) and the Dunning aug-cc-pVQZP, only for some selected structures. The T1 diagnostic parameter was also computed to check the possible multiconfigurational character (if T1 > 0.02) of the converged structures, and the expectation value of the total spin operator 〈 2 〉 was employed to estimate spin contamination. Most of the structures analyzed in this work showed multiconfigurational character and significant spin contamination, particularly at regions close to the TSs, indicating that only multiconfigurational methods such as MRCI or MS-CASPT2 provide correct wave functions. Table 9. Electronic activation energies (∆E ‡ ) and energy differences between reactants and products (∆E) for the reactivity derived from the interaction between 1 OSSO and 3 SO. All energies in kcal/mol.

Methodological details and data analyses
We use photochemical steady state calculations to estimate the abundance profiles for several sulfur species of particular importance in the Venus atmosphere. Because the photochemical lifetimes of many of the trace sulfur species are shorter than the eddy transport timescale, this approximation is valid to first order, and allows us a rapid assessment of the implications of the new chemical schemes proposed here based on our ab initio results. For non-sulfur species, long-lived sulfur species, and photodissociation rate constants, we use profiles from Pinto et al. 1  Photochemical steady state calculations are carried out from 58 to 112 km, following the temperature and total number density profiles given by Zhang et al. 27 Steady-state number densities are computed for 3 SO, 1 (SO)2, which we assume in these calculations to be primarily 1 cis-OSSO, 1 S2O, and 3 S2. A reduced set of 19 reactions involving these species is given in Supplementary Table  13. Setting production rate equal to loss rate for each of the 4 species of interest, and using the reactions in Supplementary Table 13, we arrive at the equations for steady-state number density (S1) -(S4). These equations are solved in the order presented, and the steady-state values are used as applicable. For example, the steady-state number density of 1 (SO)2, given as [(SO)2]ss and found from equation (S2), is used wherever [(SO)2] is specified in equation (S3). Loss of S2 to S4 formation is accounted for, but we are not attempting to accurately account for sulfur allotrope abundances. For this reason, S2 may be taken as a proxy for total sulfur aerosol production. 1.0 × 10 -12 molecule -1 cm 3 s -1 a Multiple scattering increases the rate constant from 5.0 x 10 -2 s -1 at 112 km to 6.2 x 10 -2 s -1 at 68 km. 27 b Updated rates of R9-R11 reactions from this work and analogous reaction of (SO)2 with H, compiled in Supplementary Table 6, do not affect the conclusions (see text in section 8.1).

S56
Equations for computing approximate steady-state number densities of SO, (SO)2, S2O and S2 given below. These expressions are approximate because i) we are not accounting for vertical transport, and ii) we are not including all possible reactions. Our objective is to get an order-of-magnitude estimate of the significance of the new branching ratios and rate coefficients determined by the ab initio calculations presented in this work, especially for reactions R3b and R13 (Supplementary Table 13). We performed several calculations to validate our approximate steady-state results versus the more complete models presented in Pinto et al. 1 Supplementary Figure 51 shows several reaction rates computed here ('calculated') compared to plots from Pinto et al. 1 We found good agreement among the rates for SO + (SO)2 and NO + (SO)2, but the reaction rate for SO + SO formation to the dimer is about a factor of 10 higher for Pinto et al. 1 compared to our calculated rate; we do not know the origin of this discrepancy. Supplementary Figure 52 shows the level of agreement for our steady-state SO and (SO)2 versus the profiles in Pinto et al. 1 Our calculated values are similar to or higher than the Pinto et al. values by at most a factor of 10. The exception, shown in Supplementary Figure 53, Figure 51. Comparison of reaction rates computed using abundance profiles from Pinto et al. 1 to plots from Pinto et al. 1 We were not able to obtain complete consistency for the 3 SO + 3 SO reaction forming the SO dimer, even with a factor of 3.3 enhancement in the low pressure rate constant to account for a CO2 atmosphere. 25 In this case both our 3 S2 and 1 S2O steady state values are low compared to Pinto et al. We do not know the cause of the discrepancy in 1 S2O. Using new UV cross sections for S2, we have increased the photodissociation rate constant for S2, which is contributing to its lower value versus Pinto et al. 1 We note that for 3 S2, and all other sulfur allotropes, condensation reactions have not been included here.

Further comments on the sensitivity of the photochemical steady state model
To evaluate the effect of updated data obtained in this work for the deoxygenation reactions of (SO)2, we run two steady state models ( Supplementary Figures 54 and 55). In the first one (Supplementary Figure 54), we use the approximate values for the NO + OSSO and O + OSSO reactions copied from the experimental value for the SO + OSSO reaction. H profile from Zhang et al. 27 was also added, including H + OSSO reaction with a high-enough rate constant of 1×10 -10 molecule -1 cm 3 s -1 . In the second one (Supplementary Figure 55), we used our updated high-level ab initio data from Supplementary Table 6. There is very little effect on the results. This is because ClS + SO and SO + OSSO completely dominate production of S2O in the present model. These results highlight the need for future ab initio and laboratory studies of reactions involving S-Cl compounds. cm 3 s -1 . For three body reactions, the low pressure limit rate constant, in units of molecule -2 cm 6 s -1 , is considered.

Origin of the rates included in the photochemical steady state model
• The photolysis rate constant for SO2 is computed at the top of the atmosphere using the solar photon flux at 1 AU from Gueymard 31 and the SO2 dissociation cross sections from the Leiden photochemical database. 32 Wavelengths shortward of 170 nm contribute a value of 2.41×10 -5 s -1 . Because we will compare our value to that of Zhang et al. 27 at 112 km, we will neglect the contribution from shortward of 170 km due to screening by the overhead column of CO2, which we estimate to be about 2×10 19 cm -2 . The B-X band from 170 to 219 nm contributes a value of 2.29×10 -4 s -1 . Scaling to Venus (a factor of 1.93) and taking the diurnal average (a factor of ½) yields the total SO2 dissociation rate constant of 2.21×10 -4 s -1 . Our value is roughly 10% larger than the value of 2.0×10 -4 s -1 obtained in Zhang et al. 27 at 112 km. We have scaled the JSO2 vertical profile from Zhang et al. 27 by the same factor at all altitudes, but this change has a minor effect.
R2 + ℎ → + • We have computed the photodissociation rate constant for SO using the dissociation cross sections from the Leiden database. 32 The band from 116 to 133 nm contributes 7.33×10 -5 s -1 , and from 190 to 235 nm contributes 4.32×10 -4 , both at 1 AU. Scaling to Venus and taking the diurnal average yields a total JSO of 4.88×10 -4 s -1 at the top of the atmosphere. Neglecting the short wavelength contribution due to overhead CO2 absorption at 112 km, we obtain 4.17×10 -4 s -1 , which is about 13% higher than the Zhang et al. 27 value of 3.7×10 -4 s -1 at 112 km.
R3a ( ) 2 + ℎ → + • The photodissociation rates for cis-and trans-OSSO were extracted from Frandsen et al. 2 . In this work, the authors calculated these photolysis rates based on theoretically predicted UV-Vis cross sections using a nuclear ensemble approach, assuming a photolysis quantum yield of 1.0, and considering a Venusian altitude of 64 km and latitude of 0º, obtaining values of 0.20 s -1 and 0.62 s -1 for cis-and trans-OSSO, respectively. In the present study, these rates were updated to 0.19 s -1 and 0.56 s -1 by considering our computed photolysis quantum yields (0.95 and 0.9 for each system, respectively).
R3b ( ) 2 + ℎ → 2 + 2 • The 1 x 10 -2 s -1 value was assumed in Pinto et al. 1 based on the assignments from the experiment of Wu et al. 33 and considering 10% of S2 production from cis-OSSO photolysis. The exact value considered in that work was 9.30 x 10 -3 s -1 . In the present study, we have considered in our analyses both the assumption from Pinto et al. (10 -2 s -1 ) and zero S2 photoproduction, being the latter based on our ab-initio results.
• The values of the rate are extracted from Zhang et al., 27 which refers to the Ph.D. thesis of Mills. 34 In Table 5A.1 on page 5.3 of said thesis, it is stated that the rate was estimated. In page 5-14 it is specified that the cross sections for the 260-340 nm range was estimated as 30 times the SO2 cross section, based on earlier work by Jones, 35 while the selected quantum yield was based on measurements by Zhang et al. 36 In Frandsen et al., 2 the authors computed the cross sections of S2O following the same approach used for cis/trans-OSSO. They also compared their spectrum with the one obtained using the 30 x SO2 cross sections estimation from F.P. Mills, being both in surprisingly good agreement. Thus, we decided to keep the values reported in Zhang et al. 27 R5 2 + ℎ → + • Using the S2 cross sections from about 230 to 290 nm in the Leiden database, 32 which are derived from Stark et al., 37 we obtained a rate constant of 7.69×10 -3 at 112 km in Venus. This is nearly a factor of 2 larger than the Zhang et al. value of 4.0×10 -3 s -1 , and results from the substantial improvement in S2 cross sections presented in Stark et al.

R6
+ → 2 + • The value of 3.5 x 10 -15 cm 3 s -1 at 298K was taken from Martinez and Herron, 38 study in which this rate was determined experimentally.
R7 + + → ( ) 2 + • Regarding the low-pressure limit of 4.4 x 10 -31 cm 6 s -1 , this value was taken from Zhang et al., 27 which refers back to the Ph.D. thesis of Mills. 34 In Table 5C.2 and page 5-34 of said thesis it is stated that the rate was taken directly from Herron and Huie. 26 In Krasnopolsky, 39 the low pressure limit rate of Herron and Huie 26 was scaled by 2.5 to account for the difference in the bath gas in their experiment (N2) compared to the Venusian atmosphere, which is mainly CO2. In Frandsen et al., 40 a factor of 3.3 was used as well to account for the different bath gas. Thus, we have adopted the latter factor of 3.3 in our simulations. • In the case of the high-pressure limit rate constant, the value of 1.7 x 10 -11 cm 3 s -1 determined using computational methods in Frandsen et al. 2 was used.
R8 ( ) 2 + → + + • For the present calculations, we neglect thermal decomposition of the SO dimer. The rate coefficient for thermal decomposition given in Zhang et al. 27 becomes smaller at higher temperature, which we consider to be problematic. We therefore set this rate constant to zero. R9 + ( ) 2 → 2 + 2 • The value of 3.3 x 10 -14 cm 3 s -1 was taken from Herron and Huie, 26 study in which experimental and modelling work was carried out to obtain the rate. In Pinto et al., 1 this rate is also included but the referenced source was Zhang et al. 27 In the latter work, Moses et al. 41 is cited as the original source, which again refers to Herron and Huie. 26 In addition, the value of the rate, depending on the study, is 3.3 x 10 -14 cm 3 s -1 or 3.0 x 10 -14 cm 3 s -1 , although 3.3 x 10 -14 cm 3 s -1 seems to be the original value in Herron and Huie. 26 In our work, an updated In addition, for reaction + 2 → + , which is only included in Supplementary Table 6, the value of 2.3 x 10 -12 cm 3 s -1 was taken from Zhang et al., 27 which refers to JPL-15 50 as the original source.
We have also considered the reaction + ( ) 2 → 2 + . A value of 3.0 x 10 -14 cm 3 s -1 can be found on Pinto et al. 1 where it was postulated that the rate for this reaction was the same as for R9. In our work, we considered both a high-enough rate constant of 1×10 -10 molecules -1 cm 3 s -1 and our updated high-level ab initio data.
Finally, we have considered the reaction 2 + ( ) 2 → 2 2 . The value of 3.3 x 10 -14 cm 3 s -1 for the rate constant was taken from Zhang et al., 27 which refers to the Ph.D. thesis of Mills 34 in which is stated that the rate was estimated.